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We investigate the superfluid properties of a Bose-Einstein condensate (BEC) trapped in a one 
dimensional periodic potential. We study, both analytically (in the tight binding limit) and numer- 
ically, the Bloch chemical potential, the Bloch energy and the Bogoliubov dispersion relation, and 
we introduce two different, density dependent, effective masses and group velocities. The Bogoli- 
ubov spectrum predicts the existence of sound waves, and the arising of energetic and dynamical 
instabilities at critical values of the BEC quasi-momentum which dramatically affect its coherence 
properties. We investigate the dependence of the dipole and Bloch oscillation frequencies in terms of 
an effective mass averaged over the density of the condensate. We illustrate our results with several 
animations obtained solving numerically the time-dependent Gross-Pitaevskii equation. 

I. INTRODUCTION 

The study of the superfluid properties of Bose-Einstein condensates (BECs) trapped in periodic potentials are 
attracting a fast growing interest. The main reason is that the control parameters of such systems are widely tunable 
in realistic experiments, allowing for the investigation of different and fundamental issues of quantum mechanics, 
ranging from quantum phase transitions [1] and atom optics [2,3] to the dynamics of Bloch and Josephson oscillations 
[4-6] . Several efforts are also focusing on the realization of new technological devices as BEC interferometers working 
at the Heisenberg limit [3] , and quantum information processors [2] . 

The dynamics of BECs in lattices is highly non-trivial, essentially because of the competition/interplay between the 
discrete translational invariance of the periodic potential and the nonlinearity arising from the interatomic interactions. 
For deep enough optical potentials, interactions induce a quantum transition from the superfluid to a Mott-insulator 
phase [1,7,8]. In this work we will study the system in a region of parameters such that its ground state stands 
deeply in the superfluid phase, with the dynamics governed by the Gross-Pitaevskii equation (GPE). Because of the 
discrete translational invariance, the excitation spectrum of the system exhibits a band structure which has several 
analogies with the electron Bloch bands in metals [9-11]. On the other hand, the coexistence of Bloch bands and 
nonlinearity allows, for instance, solitonic structures [12-14] and dynamical instabilities [15-17] which do not have an 
analog neither in metals, nor in Galilean invariant systems. 

Exact, time-dependent solutions of the GPE with an external periodic potential, Eq.(l), can be written as Bloch 
states, namely as plane waves modulated by functions having the same periodicity of the lattice. The dynamics of 
small amplitude perturbations on top of these states satisfies two coupled, linear Bogoliubov equations, which can 
be solved numerically. However, when the interwell barriers of the periodic potential are high enough, the system 
can be described in a nonlinear tight binding approximation and several important properties of the system can 
be retrieved analytically [18]. Indeed, in the nonlinear tight binding approximation the continuous Gross-Pitaesvkii 
equation can be replaced with a discrete nonlinear equation (DNL), Eq.(5), where the relevant observables of the 
system are the number of particles Ni(t) trapped in well I and the relative phases 4>i,i+i(t) = <t>i+i(t) — 4>i{t). In this 
paper, we rewrite the results derived in [18] in a more convenient form, namely in terms of two effective masses and 
group velocities. Furthermore, we compare our analytical expressions with full numerical solutions, and we extend our 
analysis to investigate the behaviour of the system at low optical potential depths, where the nonlinear tight binding 
approximation breaks down. We show that the phenomena predicted by the DNL equation (5) can be generalized to 
the case of shallow potentials, bringing new insights on the dynamics of the system. 



II. DISCRETE NONLINEAR DYNAMICS 



In the "classical" (mean field) approximation, the BEC dynamics at T = is governed by the Gross-Pitaevskii 
equation [19] 
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where g = 4:nh 2 a/m, with m the atomic mass and a the s-wave scattering length: a > (a < 0) corresponds to 
an effective interatomic repulsion (attraction). In the following we consider only a BEC with repulsive interatomic 
interactions. The external potential V ext includes the optical periodic potential Vp, which is typically superimposed 
to an harmonic (or linear) potential Vm- The periodic potential is 
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where d is the lattice spacing and 7r/d is the wavevector of the lasers in the lattice direction. The lattice spacing 
determines the Bragg momentum 
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corresponding to the boundary of the first Brillouin zone. The energy barrier between adjacent sites is expressed in 
units of the recoil energy Ep = q B /2m. From (2) we see that the minima of the laser potential are located at the 
positions xi = Id (I is an integer). Around these points, Vp « mCj 2 (x — x/) 2 /2, where Tilo x — 2^/sEp. 

The harmonic potential is Vm = m[u> 2 x 2 + io 2 y 2 + oj 2 z 2 ]/2. Since, typically, u x <C oj x , it is convenient to write the 
external potential as V ext = Vl + Vd, where the lattice potential Vl = sEp sin 2 (nx/d) + m[u> 2 y 2 + uj 2 z 2 ]/2 includes 
the transverse confining field, and the "driving" field Vd = muj 2 x 2 /2 gives the effective force acting on the center of 
mass of the condensate wave packet. 

In order to understand the basic physics of the system, we first consider the case of deep optical lattices, where 
analytic solutions can be obtained in the tight binding approximation. Then we study the behaviour of the system 
beyond the tight binding limit, solving numerically the Gross-Pitaevskii and Bogoliubov equations with arbitrarily 
shallow periodic potentials. 

As it has been previously shown in [18], when the interwell barriers are much higher than the chemical potential, 
it is possible to write the condensate wavefunction as 
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where the Wannier wavefunctions $/ are well localized in each well. The total number of atoms is Np — J2i 



J2i | 2 - Replacing this Ansatz in the Gross-Pitaevskii equation (1) and integrating out the spatial degrees of freedom, 
we find the DNL equation 



ih^ = v t + A oc 1>i-x bPi W+i + V>/-i) + c.c] fr+ 

- [K + x (I I 2 + I in+i I 2 )] in+i - [K + x (I V>H 2 + I ^-i I 2 )] Vl-i- 

The "local" chemical potential fj, l ° c is the sum of three contributions 



loc 

Vi 



J 



df 



^ mi) 2 + V L -f 2 +g\^\ 2 $f 



(5) 



(6) 



which depend on the atom number TV; explicitly through \ipi\ 2 and implicitly through the shape of the $;'s. The 
tunneling rates /Q,;±i between the adjacent sites I and I ± 1 are 
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where the on-site wavefunctions have been calculated with an average number of atoms per site, iV =| Vo | 2 ; namely 
&i(Ni) ~ $i(No) (a discussion of the validity of this approximation is in [18]). On the same line, the coefficient \ is 
given by 



X = -9 J eF$?$i±i, 

and the on-site energies V;, arising from any external potential superimposed to the optical lattice, are 
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such that Vi oc I 2 (Vi oc I) when the driving field is harmonic (linear). 

The dependence of the local chemical potential on the number of atoms is affected by the effective dimensionality 
of the condensates trapped in each well of the lattice. This can be determined comparing the interaction chemical 
potential /j,\ nt = \ipi\ g J dr <5>f and the three frequencies, cD x , oj yi u> z obtained expanding the lattice potential around 
the minima of each well Vl — m [^t{ x — x i) 2 + ^y 2 + ui 2 z 2 }/2. For instance, when u> x ,u) y ,u) z ^> n] nt , the spatial 
width of each trapped condensate does not depend (in first approximation) on the number of particles Ni in the same 
well, and the condensate wavefunction in each valley of the periodic potential is well approximated by a gaussian. 
We consider this as a OD (zero-dimensional) case: then, the nonlinear tight binding approximation (4) reduces to 
the usual tight binding approximation \&(r, t) — ^2iipi(t) [12]. The ID case arises when two frequencies are 

larger than the interaction chemical potential. In this case the system realizes an array of weakly coupled cigar- 
shaped condensates, with factorized as gaussians along the two tight directions and a Thomas-Fermi in the other 
direction. In the 2D case only one frequency is smaller than the local interaction chemical potential: we have an array 
of pancake-like condensates, with factorized as a gaussian along the tight direction and a Thomas-Fermi in the two 
other directions. The 3-D case is given by the condition [i\ nt 3> u) x ,u> y ,u> z and the wavefunction in the Zth well 4>z is 
simply given by a three-dimensional Thomas-Fermi function. The crucial point is that the effective dimensionality of 
the condensates gives a different scaling of the local interaction chemical potential (6) with the number of atoms 

= u a \ipi r, (io) 

with a = YfD> where D = 0, 1, 2, 3, \ipi\ 2 is the number of atoms in well I, and U a is a constant which does not depend 
on the number of atoms nor on the site index. When x^o "C K and D = 0, the DNL equation (5) gives the discrete 
nonlinear Schrodinger equation [12]. 



III. EXCITATION SPECTRA 



In this section we derive the Bloch and the Bogoliubov excitation spectra of the system in absence of any driving 
field (V; = 0). First we derive our results analytically in the tight binding approximation; then we solve the equations 
numerically for a wide set of parameters to extend our treatment beyond the tight binding regime. 



A. Bloch energy, Bloch chemical potential, effective masses and group velocities 

The Bloch states ^ P (x) = e lpx ^ h ^/ p (x), where ^ p {x) is periodic with period d, are exact stationary solutions of the 
Gross-Pitaevskii equation (1). Both the energy per particle e a (p) (Bloch energy) and the chemical potential /x Q (p) of 
such solutions form a band structure, so that they can be labeled by the quasi-momentum p and the band index a. 

The DNL equation (5) describes only the lowest band of the spectrum (in the following, we will consider only the 
lowest band a = 1, and we will omit, for simplicity, the band index a). Exact solutions of the DNL equation arc 
the "plane waves" ipi = V'o e l ^ pW_ ^'' i , where p is the quasi-momentum, and I is the site index (note that the ipi are 
plane waves in the discrete /-space, but do not correspond to plane waves in real space). Within the DNL equation 
framework, the energy per particle e{p) and chemical potential fi(p) corresponding to these solutions arc 

e(p) = e l ° c -2(K + 2 X N ) cos (^) ee e l ° c - cos (^) , (11) 
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where \i loc = 14°%,=^ = d(N e loc )/dN 07 with 7V = |V^o| 2 the number of atoms per well and e loc = 2U a N^ /2 /(a + 2). 
In the previous equations we have introduced the effective masses m E and m M , to emphasize the low momenta (long 
wavelength) quadratic behaviour of the Bloch energy spectrum and of the chemical potential [20] . It turns out that 
several dynamical properties of the system can be intuitively understood in terms of such effective masses. This 
approach is quite common, for instance, in the theory of metals, where = m £ . However in BEC, because of 
the nonlinearity of the Gross-Pitaevskii equation, the two relevant energies of the system, e and fi, have the same 
cos (TTp/qs) dependence on the quasi-momentum p, but different curvatures. Therefore, ^ m e , with 
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where K and x have been defined in Eqs.(7,8). Sometimes it is convenient to extend the definition of the effective 
masses to the full Brillouin zone, introducing the quasi-momentum dependent masses 
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where, following Eqs.(13,14), m £ — m £ (0) and — m„(0). 

It is also useful to introduce, with the same line of reasoning, two different group velocities, defined as 
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There is a simple, general relation between the two different group velocities (following from /i = d (N e) / dN ) 
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with, given Eqs.(13,14), > v e . The analog relation for the effective masses has been retrieved in [20]. 

Of course, the concept of effective mass, defined as the inverse of the curvature of the corresponding spectrum 
(as that of group velocity, defined as the first derivative) can be extended to shallow optical potentials, where the 
nonlinear tight binding approximation breaks down. In this case, the quasi-momentum dependence of e and pL will 
not be simply described by a cosine function, but will still remain periodic in the quasi-momentum p. In particular, 
the value p where m £ (p) changes sign (corresponding to d 2 e/dp 2 = 0) will be greater than gs/2 and will in general 
not coincide with the momentum where m f _ l (p) changes sign (corresponding to d 2 /j,/dp 2 = 0). 

The presence of the two different effective masses (group velocities) raises an important problem: which effective 
mass (group velocity), and how, enters in the dynamical properties of the system? For instance, we anticipate that 
the current carried by a Bloch waves with quasi-momentum p is po v s (p), where po is the average particle density; m M , 
on the other hand, plays a crucial role in the Bogoliubov spectrum. To conclude this subsection, we notice that the 
Bloch states are not the only stationary solutions of the Gross-Pitaevskii equation. Because of nonlinearity, indeed, 
periodic solitonic solutions can also appear for a weak enough periodic potential, introducing new branches in the 
excitation spectra [21,22]. 



B. Bogoliubov dispersion relation 



In this subsection we study the Bogoliubov spectrum of elementary excitations. This describes the energy of 
small perturbations with quasi-momentum q on top of a macroscopically populated state with quasi-momentum p 
(stationary solution of Eq.(l)). To be explicit, let us consider first the case in which the radial degrees of freedom y, z 
are integrated out. The wavefunction along the x direction can be written as 
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Because of the periodicity, the Bogoliubov amplitudes can be written as Bloch waves [i.e., {u, v} pq (x) = 
exp(iqx/h){u,v} pq (x)], where q is the quasi-momentum of the excitation and {u, v} pq (x) are periodic functions. 
The subscript {pq} indicates that both the amplitudes u, v and the excitation frequencies uj pq depend on the quasi- 
momentum p of the carrying wave and on the quasi-momentum q of the excitation. 
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In terms of the periodic functions ^, u and v, the Bogoliubov equations take the form 
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where n is the 3D-average density and f_ d ^ 2 \^ p \ 2 dx — 1. Equations (21,22) can be solved numerically in a very 
efficient way working with the Fourier components of it and v. 



(a) s=20, gn=0.5E R , p=0.4q B (b) s=20, gn=0.5E R , p=0.7q t 




q/q B q/q B 

FIG. 1. Numerical solutions of Eqs.(21,22) for s = 20 and gn = 0.5Er (green dots); analytic solution Eq.(23) in the tight 
binding approximation (solid blue line); analytic solution of Eq.(23) where m M is replaced with m e (dashed red line). The 
quasi-momentum of the carrying wave is p — OAqs in (a) and p — 0.7gs in (b). 



In the tight binding regime, the Bogoliubov analysis corresponds to perturbing the large amplitude wave as ipi — 
[ipo + e*( pW ~ Alt '/ ?i , with (5-0; = 'J2^U q e^ qld ^ h ~ UJpqt \ Retaining only first order terms with respect to Sip, we get two 
coupled linear equations analogous to (21) whose eigenvalues can be calculated analytically [18]. The general solution 
(for any effective dimensionality of the system: D = 0, 1, 2, 3) is 
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with the chemical potential given by fJ,(p,N ) = [i loc - cos (^J (see Eq.(12)), and [i\ oc = U a | ipi \ 4 /<- 2+D \ 

For D = and in the limit \ — 0; we recover the well-known results for the discrete nonlinear Schrodinger equation 
[16,24]. Equation (23) has been first written in [18] in terms of the parameter of the DNL equation Eq.(5), while, for 
small q and arbitrary p, has been derived in [17] for arbitrary values of s. 

In Fig. 1, we compare the analytic results (dots) with the numerical solution of Eqs.(21,22) (solid line), for a 
system in tight binding regime. In the numerical calculations, the effective masses have been obtained from the 
curvatures of the Bloch energy and chemical potential spectra, while the term No-Jjj^ has been evaluated from the 
density dependence of the chemical potential. As it has already been noted in [20], effects related with the difference 
between the two effective masses in the Bogoliubov spectrum of a condensate at rest (p = 0) are usually negligible. 
On the contrary, such difference becomes important when the condensate moves with a large quasi-momentum, as 
shown in Fig. 1(b). 
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IV. SOUND WAVES & INSTABILITIES 



The small q (large wavelength) limit of the Bogoliubov dispersion relation becomes 
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(we assume, for the moment, that -^--^-Nqcos^p/ q B ) > 0). The linear behaviour in q indicates that the system 
supports (low amplitude) sound waves, propagating on top of the large amplitude traveling wave Hf p with velocity 
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where the "chemical potential group velocity" w M has been defined in Eq.(18), and the "relative sound velocity" c is 
defined as 
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The two different velocities v s> ± correspond, respectively, to a sound wave propagating in the same and in the opposite 
direction of the large amplitude traveling wave. As we have already noticed, and we will discuss again in the next 
section, is different from (it is larger than) the actual velocity of the large amplitude wave, see Eq.(19). 

We remark that, contrary to the case of a Galilean invariant system (s = 0), the sound velocity depends on the 
quasi-momentum p. Moreover, v s depends on the effective dimensionality of the condensates, since (from Eq.(10,12)) 
-J^-Nq ~ a U a Nq^ 2 . In the limit a = 2, p — > and m e , — > m we get the sound velocity in the uniform case. 
° J Th e system is energetically unstable if there exist any aj pq < 0. In the limit s = 0, this corresponds to a group velocity 
larger than the sound velocity (Landau criterion for superfluidity). When the system has a discrete translational 
invariance (s > 0) the condition for this instability is obtained from the Bogoliubov excitation spectrum Eq.(23). 
Then, we have that the system is not superfluid when u> pq < 0, corresponding to 
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This result should be compared with the well known Landau criteria for an homogeneous system (s = 0), stating that 

u = §| = §p being the group velocity of the condensate, and 



the superfluid is energetically unstable when v 2 > c 2 , 



— S&N the sound velocity 

There is a further dynamical (modulational) instability mechanism associated with the appearance of an imagi- 
nary component in the Bogoliubov frequencies, which disappears in absence of interatomic interactions, or in the 
translational invariant limit (if a > 0). The onset of this instability in the tight binding regime, coincides with the 
condition 
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The dynamical instability drives an exponentially fast increase of the amplitude of the (initially small) fluctuations of 
the condensate. Since the initial phases and amplitudes of the fluctuation modes are essentially random, their growth 
induce a strong dephasing of the condensate, and dissipates its translational kinetic energy (which is transformed 
in incoherent collective and single particles excitations). The unstable modes q grow with a time scale given by the 
imaginary part of the excitation frequency 
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We remark here the different scaling of the energetic and dynamical instability with the interatomic interactions. 
Decreasing the scattering length, the sound velocity decreases, and smaller and smaller group velocities can breakdown 
the superfluidity of the system (in the limit a = 0, the sound velocity c = 0: the non interacting condensate is 
always energetically unstable for an arbitrary small group velocity). On the other hand, the dynamical modulational 
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instability criteria does not depend on the scattering length. This apparent paradox is simply solved noticing that 
the growth time of the unstable modes, Eq.(29), actually depends on interactions, and diverge when the scattering 
length vanishes (r — > oo when a — > 0). Therefore, a noninteracting condensate is always dynamically stable. There is 
a further point to remark: if we consider a condensate moving with an increasing velocity, the system always becomes 
first energetically unstable, then it hits the dynamical instability. As a matter of fact, however, in real experiments 
the energetic instability can grow quite slowly (and at zero temperature only in presence of impurities [15]), so that 
the dominant dephasing mechanism is given by the modulational instability. This aspect can be highlighted also 
with numerical experiments, studying, for instance, Bloch oscillations of a condensate with the interactions switched 
off. In this case, even though the system is energetically unstable, it remains coherent over many oscillations. If the 
interatomic interactions are switched on, however, the system dephase rather quickly, the dephasing occurring when 
the quasi-momcntum of the condensate is in the dynamical unstable region of the Bloch spectrum. We have done 
such numerical experiment, and results are shown in Movies 4 and 5. Of course, our prediction can be tested in real 
experiments, tuning the scattering length with a Feshbach resonance. 

To summarize, the tight binding approximation predicts the arising of the dynamical instability (complex excitation 
frequency) for p = qB^- We point out that p = qs/2 also corresponds to the quasi-momentum where, in the tight 
binding regime, the effective masses m e (p) and m^{p) change sign. A system with a negative effective mass and 
positive scattering length can be, roughly speaking, seen as equivalent to a system with a negative scattering length 
and positive effective mass. It is well known that a BEC having a negative scattering length is dynamically unstable, 
and, therefore, such parallelism could be proposed to give a simple explanation of the instability. However, we will 
see that this coincidence between the arising of dynamical instabilities and the inversion of sign of the effective mass 
does not take place at lower optical potentials (see Fig. 3). 

Let us concentrate now on the behaviour of the excitation frequencies for shallow optical potentials, where the tight 
binding expression derived in (23) is not applicable. For small optical potential depths, the Bogoliubov equations 
have to be solved numerically and the results show a more complicated behaviour. In Movies 2(a-c), we show the 
numerical solutions of the ID Bogoliubov equations for three different values of s (s = 1,2 and 5). In those movies, 
we plot the real and imaginary part of uj pq as a function of q and we vary p in time. 
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FIG. 2. MOVIES: Real and imaginary part of u) pq as a function of q for different values p, for gn — 0.5Er and s = 1 (a), 
s = 2 (b) and s — 5 (c). 



We point out a series of differences with respect to the tight binding regime: 

- the complex frequencies appear at the boundary of the first Brillouin zone (q = qs) for a value of p > qs/'i (dots 
in Fig. 3) and they reach the center of the zone (q = 0) for a higher value of p (orange region in Fig. 3). 

- the range of momenta p where the frequencies are complex for some g, but real around q = 0, decrease by increasing 
s. In the tight binding limit this range vanishes; 

- in the limit of our numerical accuracy, which is due to the discrete sampling of p and q, we found that the value of 
p where the effective mass m e changes sign (squares in Fig. 3) corresponds to the value of p at which the frequencies 
with non vanishing imaginary part reach q w 0. In the tight binding approximation, this appears explicitly through 
the term cos(Trp/qB)/m £ under the square root. 

We would like to remark two important results arising from our study of the excitation spectra. First, as shown 
in Fig. 3, we found the onset of the dynamical instability for values of the quasi-momentum where the effective mass 
m e (p) is still positive. Second, the range of momenta where the system has a positive effective mass and, at the same 
time, is dynamically unstable, increases by decreasing s (keeping in mind that the amplitude of the imaginary part 
vanishes for s — ► or gn — > 0, which implies that the growth in time of the instability diverges both for uniform 
interacting systems and for an ideal gases in optical lattices). So, one can study the behaviour of the system at low 
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s to distinguish between two possible dephasing mechanisms, one due to the sign of m e , the other one due to the 
dynamical instability, as it will be extensively explained in Sect. VI. 

We conclude this section remarking that various important aspects of the physics of energetic and dynamical 
instabilities of a BEC in a periodic potential have been studied in [12,15-17]. 
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FIG. 3. Results from the numerical solution of Eqs. (21,22) for gn — 0.5Er as a function of s. Dots: value of 
the quasi-momentum p where the excitation frequencies uj pq (q = ±<?b) become complex; orange region: value of the 
quasi-momentum p where complex frequencies are found around q — 0; squares: value of the quasi-momentum p where 
the effective mass m e changes sign. The dotted lines are a guide to the eye. 

V. NEWTONIAN DYNAMICS 

Using the results in [18], we can now rewrite the dynamics of a BEC wave-packet in terms of the energy effective 
mass. For the BEC wave-packet we use the following ansatz 



(30) 



where £(t) and a{t) are, respectively, the center and the width (in lattice units) of the wavepacket, P{t) and 8{t) 
their associated momenta and JC{a) a normalization factor such that J^i^i = Nt (with \ipi\ 2 = N{). The function 
/ is generic, for instance, we can choose f{X) = e- x ' 2 or f(X) = (1 - X 2 ) 1 /" (with -1 < X < 1) to describe, 
respectively, the dynamics of a gaussian or a Thomas-Fermi wavepacket. The equations of motion of the collective 
variables £(t), a(t),p(t), S(t) have been obtained in [12,18]. With Vj = Vtj 2 (Q = md 2 uj 2 /2), and neglecting the 
dynamics wave-packet width dynamics (<r(t) = 0), we hnd that the group velocity £ and the effective force acting on 
the center of mass of the wavepacket are given by 



7r^ \ m F 



hP=- 



dVd 
9£ 



(31) 
(32) 



where V d = fl(£_ 2 + cr 2 ^) with X x = J dXf 2 (X) and 1 2 = J dXX 2 f 2 (X). Since the effective masses depend on 
the local (on-site) density, we have to introduce an effective mass averaged over the local density of the condensate 
wavepacket 



1 



(33) 



with, according to Eq.(13), m~ 1 (Ni) = (27r 2 /g^) (/<" + 2\Ni). We summarize here the most important results, written 
in term of the effective mass m e : 

(i) in the case of an homogeneous system (Vd = 0, Ni = const.), the tunneling rate is given by 



1 
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(34) 
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(ii) the frequency of small amplitude oscillations of the wavepacket driven by an harmonic field V; oc I 2 i 



is 



'x 




) 



(35) 



(Hi) if the driving field is linear Vi = mGdl, we have simple Bloch oscillations with 




irm G 



t 



) 



(36) 



This analysis does not take into account possible dcphasing mechanisms as those investigated in the previous section. 
In the collective coordinate approach, such dephasing mechanisms can be described including the dynamics of the 
width of the wavepacket a(t) and of the corresponding momentum S(t) [12]. 

VI. NUMERICAL EXPERIMENTS ON BLOCH OSCILLATIONS, DIPOLE OSCILLATIONS AND FREE 

EXPANSIONS IN THE LATTICE 

In this section, we discuss some numerical simulations of the Gross-Pitaevskii equation in order to illustrate the 
phenomena described in the previous sections. We first consider Bloch oscillations (Sect. VI A): we create the con- 
densate in a harmonic trap superimposed to the lattice and then switch off the harmonic trap and replace it with a 
linear potential. We expect that the BEC oscillates periodically in space (Bloch oscillations). 

The second numerical simulation consists in creating the condensate in a harmonic trap superimposed to the lattice, 
and suddenly displacing the center of the harmonic trap (Sect.VIB). This experiment has already been studied 
theoretically [16] and performed experimentally in [6,23]. We discuss it again, generalizing the previous results to the 
case of shallow optical lattices. 

The third numerical simulation consists in creating a condensate in a harmonic trap superimposed to the optical 
lattice and, thereafter, switching off the harmonic trap in the lattice direction, letting the condensate expand in the 
periodic potential (Sect. VIC): for values of the the mean-field energy large enough (with a fixed height s of the 
interwell energy barriers), the wave-packet is self-trapped [12] and the spreading of the wave packet does not occur. 

In all cases, for an interacting BEC, we have found some sort of self-trapping and dynamical instabilities for some 
values of the depth of the periodic potential or the initial conditions of the BEC wave packet. For instance, in 
the dipole oscillation experiment, the condensate may stop on one side of the harmonic potential being unable to 
complete the oscillation. It is useful to look at the dynamical evolution of the relative phases of condensates trapped 
in neighboring wells and the BEC evolution in momentum space. There is a clear correspondence between the 
distribution in momentum space and that in quasi-momentum space: the quasi-momentum distribution of the Bloch 
state typ(x) = e lpx / h 'fy p (x) with quasi-momentum p is 5(p); its momentum distribution is | ^2 e cgS(p + 2£qB)\ 2 , where 
I are integers and where the q are the Fourier coefficient of the periodic function $ p (x). An analogous relation is also 
valid when the condensate is not in a well-defined Bloch state, but in a superposition of Bloch states of the first band. 
In this case, the width of the peaks of the momentum distribution will be equal to the width of the quasi-momentum 
distribution. In the following we will work with the momentum distribution, which is simply obtained from the Fourier 
transform of the condensate wavefunction. 



Bloch oscillations can be explained in very simple terms. In the presence of a linear potential superimposed to the 
optical lattice, the behaviour of a particle with quasi-momentum p is described by the equation of motion p(t) = Ft, 
where F is the constant force due to the linear potential. Since the velocity is given by v g = de(p)/dp, when the 
effective mass is negative, the particle will respond to a positive (negative) force with a negative (positive) acceleration. 
Since the energy band e(p) is periodic in p, this will result in periodic oscillations in coordinate and velocity space. 

This simple explanation, even neglecting important effects like Landau-Zener tunneling to a higher band, provides 
a useful model to interpret experiments with electrons [25] , with cold atoms [26] and with Bose-Einstein condensates 
[4,5]. However, if interactions in the condensate play a mayor role, the scenario can change dramatically. First of all, 
the momentum distribution ^>(p, t) will not evolve just like \&(p(i)), as it approximately happens for non interacting 
systems, but will also spread. Furthermore, in the presence of interactions, it might happen that the condensate 
gets dephased and, after a short while, the oscillations stop (see Movie 4). For the situation described in the movie 



A. Bloch oscillations 
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(s = 10, for which the tight binding approximation works well), the dephasing process begins when the center of the 
momentum distribution reaches p = <7b/2. This point corresponds both to the on-set of the dynamical instabilities 
and the inversion of the effective mass. Since the momentum distribution has a certain width, one could think that 
the oscillations stop because the sign of the effective mass is not the same for the whole condensate. An alternative 
interpretation relies on the onset of dynamical instabilities. 
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FIG. 4. MOVIE: Bloch oscillations for s = 10 and gn — 0.5Er. Time evolution of the spatial density (upper plot), of the 
relative phases (middle plot) and of the momentum distribution (lower plot). 
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FIG. 5. MOVIE: Bloch oscillations with s = 10, gn = (green line) and gn = 0.5Er (blue line). Time evolution of the 
spatial density (upper plot), of the center of mass (middle plot) and of the momentum distribution (lower plot). The red dotted 
line indicates the onset of dynamical instabilities in the interacting case. 

In order to highlight the correct interpretation, we study the Bloch dynamics of a non interacting condensate, which 
is always dynamically stable. The initial spatial width is chosen in order to get about the same momentum distribution 
as in the interacting case, in order to have similar effective mass effects. More specifically, since in the interacting case 
the width of the momentum distribution increases slowly, while in the non interacting case it remains almost constant, 
we choose the initial conditions so that the two momentum distributions will be similar at the "critical point" , where 
(P) - Qb/2. 
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The direct comparison is shown in Movie 5. We observe, in the non-interacting case, regular, perfectly periodic 
Bloch oscillations, in spite of the finite width of the momentum distribution. This clearly shows that, in the interacting 
case, the onset of decoherence is due to the dynamical instability. 



B. Dipole oscillations 



Dipole oscillations consist in the motion of the condensate at the bottom of the harmonic trap. The average velocity 
is periodic in time and the momentum distribution, showing the characteristic peaks due to the optical lattice, also 
oscillates periodically in time at the bottom of the band. During the time evolution, the phase differences between 
neighbouring condensates remain locked over the whole condensate. 

For a given set of parameters corresponding to small displacements, small interactions or small optical potential 
depths, dipole oscillations remain periodic, with the condensate locked in phase. Instead, increasing one of these 
quantities, we find that the oscillations get dephased during the time evolution, or even stop before the condensate 
reaches the bottom of the harmonic potential. For the seek of comparison we display in Movie 6(a-c) the evolution of 
the density, of the phase difference and of the momentum distribution for the following sets of parameters. For fixed 
interactions (gn = 0.01-Er) and harmonic trap (hcu x = 0.004_Er), we choose: 

(a) s = 3, xo = 3d: oscillations; 

(b) s = 3, Xq = 9d: broken oscillations; 

(c) s = 10, Xq = 9d: broken oscillations. 
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FIG. 6. MOVIES: Dipole oscillations for gn = OMEr, s = 3, x = 3d (a); gn = 0ME R , s = 3, x = 9d (b); gn = O.OIEr, 
s = 10, xo — 9d (c); gn — 0.5Er, s = 1, Xo = 9d (d). Time evolution of the spatial density (upper plot), of the relative 
phases (middle plot) and of the momentum distribution (lower plot). In movie (d), we indicate with a red dotted line the 
quasi-momentum where the dynamical instabilities arise and with a green dotted line the quasi-momentum where the effective 
mass changes sign. 



Looking at the phase difference between neighbouring condensates, we find that when the condensate oscillation is 
interrupted, the phases get scrambled. This corresponds to a randomized flux of atoms which arc not able anymore to 
flow coherently downwards the potential. The evolution of the momentum distribution suggests that this phenomenon 
happens when the condensate reaches the instability region, given in the specific cases by p greater than 

To further explore this interpretation, we choose a shallow optical potential such that there is a broad range of 
p where the effective mass m e (p) is positive and at the same time the system is dynamically unstable (see Fig. 
3). We increase the nonlinear interaction parameter to get a relevant imaginary part of the excitation frequencies, 
otherwise the time scale where instabilities manifest themselves is too long. In Movie 6(d) (lower plot), we indicate 
with a red dotted line the quasi-momentum where the dynamical instabilities arise and with a green dotted line 
the quasi-momentum where the effective mass changes sign. We actually observe the first signatures of decoherence 
when the momentum distribution is contained between the two lines, indicating than the decoherence happens in 
correspondence of the dynamical instability point. We conclude this section mentioning that experimental evidences 
of dynamical instabilities are reported in [23] . 
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C. Expansion in the lattice 



After creating the condensate in the harmonic trap superimposed to the lattice, we switch off the harmonic trap and 
let the condensate, which is initially at rest, expand. During the expansion, the current of atoms is from the inside 
to the outside of the cloud and the phase differences increase, being positive for x < and negative for x > 0. In [12] 
the occurrence of self-trapping has been predicted in the tight-binding: when interactions are larger than a critical 
value, the width of the wave packet does not continue to increase with time (as for vanishing or small interactions) 
and the wave packet remains localized around a few sites. A similar nonlinear self-trapping occurs in the two-site 
problem [27]. 
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FIG. 7. MOVIES: Expansion of a BEC wave-packet for gn = 0, s = 5 (a); gn = 0.01-Er, s = 5 (b); Time evolution of the 
spatial density (upper plot), of the relative phases (middle plot) and of the momentum distribution (lower plot). 



Increasing the interactions, the system enters in the self-trapped regime as shown in Movie 7(b). If interactions 
are strong enough, we see that after a first stage (whose duration depends on the strength of the interactions) the 
expansion stops and the condensate evolves as a random flow of atoms between the condensates localized at the 
bottom of the different potential wells, indicating the onset of a new dynamical instability. In this case, however, 
a Bogoliubov-like stability analysis is much more problematic because of the non trivial temporal evolution of the 
condensate wave-packet. A possible, though approximate, approach is to write down an effective Hamiltonian of the 
system in terms, for instance, of the collective coordinates introduced in Sect.V. Such Hamiltonian would contain a 
limited number of degree of freedom, making the stability analysis a much easier task. This is the approach followed 
in [12] to study the dynamics of an expanding condensate in the discrete nonlinear Schrodinger equation framework 
(with rrifi — m £ ). Within this approach one recovers, in a unified framework, the critical values of the parameters 
for the self-trapping conditions of a wavepacket of finite width initially at rest, and the onset of the modulational 
instability of a Bloch wave discussed in the previous sections. For instance, considering a gaussian wavepacket with 
initial width cto and quasi- momentum the collective coordinates approach predicts the onset of self-trapping at a 
critical value of the parameter A = U^NtJ^K [12]. When cos (po) > 0, the critical value is A c w 2y / 7rero cos (po)', when 
cos (po) < 0, the critical value is A c 2y / 7r | cos (po) | /(Jq. We remark that when the width of the wavepacket is very 
large (ao — * oo), A c — > oo if cos (po) > (and the system is always dynamically stable), while A c — > if cos (po) < 
(and the system is always dynamically unstable), recovering the findings of Sect. IV. 

The study of the dynamical instabilities of a condensate trapped in an periodic potential is quite a rich problem, 
and deserves further investigations. As we have mentioned, this is connected to the general problem of the stability 
of a non stationary state, which includes for instance also the propagation of sound waves in the non linear regime. 
First experimental results on the self-trapping with weakly coupled BECs are reported in [28,29]. 
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VII. CONCLUSIONS 



The Gross-Pitaevskii dynamics of a Bose-Einstein condensate trapped in a deep periodic potential can be studied 
in terms of a discrete, nonlinear equation. This mapping allows a clear and intuitive picture of the main dynamical 
properties of the system, which can be calculated analytically. We have calculated the effective masses of the system, 
connected to the Bloch energy and chemical potential spectra. We have calculated the Bogoliubov dispersion relation, 
and studied the sound velocity and the appearance of energetic and dynamical instabilities. We have generalized these 
concepts to the case of shallow optical lattice, which requires a numerical solution and provides complementary insight 
in the understanding of the problem. Both in the tight binding limit and in the case of shallow optical potential, we 
have investigated in detail the arising of dynamical instabilities, which seem to be the main mechanism of dcphasing 
of the condensate in Bloch oscillation and dipole oscillations experiments. 

Acknowledgements. We thank M. Kramer, L. Pitaevskii and S. Stringari for interesting discussions. This work has 
been partially supported by the DOE. 

Note added in Proofs: An equation similar to the DNL equation (5) has been derived in M. Oster, M. Johansson, 
and A. Eriksson, Phys. Rev. E 67, 056606 (2003), to describe the dynamics of an electric field in an array of coupled 
optical waveguides embedded in a material with Kerr nonlincarities. 
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